program main
  use iso_c_binding, only: c_ptr
  implicit none
  include "uhm/wrapper/fort/parameters.f90"

  character(len=256) :: argv
  character(len=32) :: &
       opt_blk ='-blk', text_blk ='setting blocksize', &
       opt_rhs ='-rhs', text_rhs ='setting # rhs', &
       opt_time ='-time', text_time ='setting time measuring level'

  integer :: iel, rhs, blk, time, i, j, k, istat

  integer, parameter :: ids = 1
  integer, parameter :: n_nodes = 7, n_elts = 2
  integer, dimension(n_nodes, n_elts), parameter :: &
       nodes = reshape((/&
       3,4,5,  7, 9,8,  1,   &
       4,6,5, 10,11,9,  2/), (/7,2/))

  integer, parameter :: dof = 15
  integer, dimension(n_nodes), parameter :: &
       dofs    = (/1,1,1,  3,3,3,  3/), &
       weights = (/0,0,0,  0,0,0,  0/), &
       kinds   = (/0,0,0,  0,0,0,  0/)

  real*8, dimension(:,:), allocatable :: A, B, X
  real*8 :: residual

  logical :: iverbose

  type(c_ptr), dimension(2) :: solver

  call get_command(argv)

  call uhm_initialize(argv)

  ! Fortran send string in a weird way
  call uhm_option_begin
  call uhm_get_option_int(trim(opt_blk),  trim(text_blk),  256, blk)
  call uhm_get_option_int(trim(opt_rhs),  trim(text_rhs),  1, rhs)
  call uhm_get_option_int(trim(opt_time), trim(text_time), UHM_TIME_LEVEL_FRONT, time)
  call uhm_option_end

  call uhm_get_verbose(iverbose)
  call uhm_time_begin(time)

  call uhm_direct_lu_nopiv_initialize( &
       UHM_DOUBLE, rhs, blk, solver)

  allocate(A(dof,dof), B(dof,rhs), stat=istat)
  if (istat /= 0) stop -1

  call uhm_solver_clear(solver)

  call uhm_solver_copy_in_begin(solver)  
  do iel=1,2
     call srand(100)

     do j=1,dof
        do i=1,dof
           A(i,j) = rand()
           if (iverbose) then
              write(*,*) 'A(',i,',',j,') = ', A(i,j) 
           end if
        end do
     end do
     do j=1,rhs
        do i=1,dof
           B(i,j) = rand()
           if (iverbose) then
              write(*,*) 'B(',i,',',j,') = ', B(i,j) 
           end if
        end do
     end do
     if (iverbose) then
        write(*,*) 'iel = ', iel, 'nodes = ', nodes(:,iel)
     end if
     call uhm_solver_copy_in_ab( solver, &
          iel, ids, &
          nodes(:,iel), ids, dofs, weights, kinds, n_nodes, &
          dof, A, B)
  end do
  call uhm_solver_copy_in_end(solver)

  deallocate(A, B, stat=istat)
  if (istat /= 0) stop -1

  call uhm_solver_create(solver)

  call uhm_direct_solver_decompose(solver)
  call uhm_direct_solver_solve(solver)
  call uhm_direct_solver_check(solver, residual)

  if (iverbose) then
     call uhm_solver_report(solver)
  end if

  write(*,*) 'residual = ', residual

  allocate(X(dof,rhs), stat=istat)
  if (istat /= 0) stop -1

  do iel=1,2
     X = 0.d0
     call uhm_solver_copy_out_el( solver, &
          iel, ids, dof, X);
     if (iverbose) then
        write(*,*) 'iel = ', iel
        do i=1,dof
           write(*,*) 'i= ', i, 'X = ', X(i,:)
        end do
     end if
  end do

  deallocate(X, stat=istat)
  if (istat /= 0) stop -1

  call uhm_solver_clear(solver)
  call uhm_solver_flush(solver)

  call uhm_solver_finalize(solver)

  call uhm_time_end
  call uhm_finalize

end program main
